lcg.sample <- function(nsamples=5000, r=10) {
  px <- c(0, cumsum(lcg.cutdist(lcg.dist.x, r, 1:r)))
  v <- findInterval(runif(1:nsamples), px)
  return(v)
}

lcg.sample1 <- function(nsamples=5000, r=10) {  
  trans <- function(r) {
    blue<-c(1)
    if (r > 1) {
      for (v in 2:r) {
        if (sample(1:v,1) == 1)
          blue <- c(blue, v)
      }
    }
    return(c(sample(blue,1), length(blue)))
  }
  
  transmatrix <- data.frame(seller=numeric(nsamples),
                            clique.size=as.numeric(nsamples))  
  for (i in 1:nsamples) { transmatrix[i,] = trans(r) }
  return(transmatrix)
}

lcg.cutdist <- function(distarray, n=1, x=1) {
  pos <- n*(n-1)/2
  probs <- distarray[(pos + 1):(pos + n)]
  return(probs[x])
}

lcg.dist.gen1r <- function(n, debug=FALSE) {
  v <- 1:n
  p <- -log(1:n)
  q <- log(1 - (1:n)**-1)
  d <- rep(0, n*(n+1)/2)
  pos <- v*(v-1)/2

  ss <- 0
  nn <- 2**(n-1)
  bt <- c(1, rep(0,n-1))
  btc <- rep(1,n)
  for (k in 0:(nn - 1)) {
    s <- k # 某组合的整数表示
    i <- 2 # 组合整数的第i个二进制数位
    t <- 0
    h <- 1
    while (i <= n) {
      if (s %% 2 > 0 ) {
        t <- t + p[i]
        h <- h + 1
        bt[i] <- 1  # 选中
      } else {
        t <- t + q[i]
        bt[i] <- 0  # 没中
      }
      btc[i] <- h
      s <- trunc(s/2); # 截掉最后一个二进制数位
      i <- i + 1
    }
    t <- exp(t)

    for (v in 1:n) {
      for(i in 1:v) {
        if (bt[i] > 0) {
          d[pos[v] + i] <- d[pos[v] + i] + t/btc[v]
        }
      }
    }
  }
  return(d)
}


lcg.dist.gen2r <- function(r) {
  #使用迭代计算生成LCG分布
  px <- rep(as.numeric(0), r*(r+1)/2)
  pm <- rep(as.numeric(0), r*(r+1)/2)
  d <- matrix(nrow=r, ncol=r)
  d[1,1] <- 1
  h <- 1
  px[h] <- 1
  pm[h] <- 1
  for (n in 2:r) {
    for (x in 1:n)
      d[n,x] <- d[n-1,1]*(n-1)/n/n

    if (n > 2) {
      for (m in (n-1):2) {
        for(x in (n-1):1) {
          d[m,x] <- d[m,x]*(n-1)/n + d[m-1,x]*(m-1)/m/n
        }
      }
    }

    d[1,1] <- d[1,1]*(n-1)/n
    d[1,2:n] <- 0

    if (n >2) {
      for (m in 2:(n-1)) {
        d[m,n] <- (m-1)*d[m,1] - sum(d[m,2:(n-1)])
      }
    }
    for (i in 1:n) {
      h <- h + 1
      px[h] <- sum(d[1:n,i])
      pm[h] <- sum(d[i,1:n])
    }
  }
  return(list(x=px,m=pm))
}


lcg.c.compile <- function() {
writeLines("
#include <math.h>
#include <string.h>
#include <stdlib.h>

// Search all space with combination, O(2^n)
void lcg_gen_dist1(int *n, double* d) {
   if (*n <=0 || *n > 48)
    return;
  double *p, *q;
  int i,j;
  double t;
  p = malloc(*n*sizeof(double));
  if (p == NULL)
    goto END;
  q = malloc(*n*sizeof(double));
  if (q == NULL)
    goto END;  
  *p = 0; *q = 0;
  for (i = 2; i <= *n; i++) {
    t = -log(i);
    p[i-1] = t;
    q[i-1] = log(i - 1) + t;
  }
  memset(d, 0, *n*sizeof(double));  
  unsigned int nn;
  unsigned long long  k;
  nn = 1 << *n;
  for (k=1; k <= nn; k+=2) {
    unsigned long long s;
    unsigned int h;
    s = k>>1;
    h = 1; 
    t = 0;
    i = 1; 
    while(i < *n) {
      if (s & 0x1) {
	t += p[i];
	h+=1;
      } else {
	t += q[i];
      }
      i++;
      s >>= 1;
    }
    t = exp(t);
    double h1 = 1.0/h;
    for (i=0; i < *n; i++) {
      s = k;
      for(i = 0; i < *n; i++) {
	if (s & 0x1) {
	  d[i] += t*h1;
	} 
	s >>= 1;
      }
    }
  }
 END:
  if (p) 
    free(p);
  if (q) 
    free(q);  
} 
// Compute distribution with iterated formula
typedef long double numeric;
void lcg_gen_dist2(int *pr, double *px) {

  int r;
  int m,n,x;
  numeric *d, *p0, *p1, *p1n, *q,*qn, t;

  r = *pr;

  d = malloc(sizeof(numeric)*r*r);
  if(d == NULL)
    goto ERROR;
  memset(d, 0, sizeof(numeric)*r*r);

  p1 = d;
  *p1++ = 1;

  *px++ = 1;

  for (n=2; n<=r; n++) {
    numeric n1n = ((numeric)(n-1))/n;

    p1 = d + r*n - r;
    p1n = p1 + n;
    t = *(p1-r)*n1n/n;
    while (p1 < p1n) {
      *p1++ = t;
    }

    if (n > 2) {
      for (m= n-1; m > 1; m--) {
	p1 = d + r*m - r;
	p1n = p1 + n - 1;
	p0 = p1 - r;
	numeric m1m = ((numeric)(m-1))/m/n;
	while (p1 < p1n) {
	  *p1 = *p1*n1n + *p0*m1m;
	  p1++; p0++;
	}
      }
    }
    p0 = d;
    *p0 = *p0*n1n;

    if (n > 2) {
      p1 = d + r;
      p1n = d + r*(n-1);
      m = 2;
      while (p1 < p1n) {
	t=0;
	q = p1 + 1; 
	qn = q + (n-1);
	while (q < qn)
	  t += *q++;

	*(p1 + n-1) = *p1*(m-1) - t;
	p1 += r;
	m++;
      }
    }

    p1 = d;
    p1n = d + n;
    qn = d + r*n;
    while (p1 < p1n) {
      q = p1;
      t = 0;
      while(q < qn) {
	t += *q; 
	q +=r;
      }
      *px = t;
      p1++;
      px++;
    }
  }  

  if (d) 
    free(d);
 ERROR:
  return;
}

void lcg_gen_dist_m(int *pr, double *px) {

  int r;
  int m,n,x;
  numeric *d, *p0, *p1, *p1n, *q,*qn, t;

  r = *pr;

  d = malloc(sizeof(numeric)*r*r);
  if(d == NULL)
    goto ERROR;
  memset(d, 0, sizeof(numeric)*r*r);

  p1 = d;
  *p1++ = 1;

  *px++ = 1;

  for (n=2; n<=r; n++) {
    numeric n1n = ((numeric)(n-1))/n;

    p1 = d + r*n - r;
    p1n = p1 + n;
    t = *(p1-r)*n1n/n;
    while (p1 < p1n) {
      *p1++ = t;
    }

    if (n > 2) {
      for (m= n-1; m > 1; m--) {
	p1 = d + r*m - r;
	p1n = p1 + n - 1;
	p0 = p1 - r;
	numeric m1m = ((numeric)(m-1))/m/n;
	while (p1 < p1n) {
	  *p1 = *p1*n1n + *p0*m1m;
	  p1++; p0++;
	}
      }
    }
    p0 = d;
    *p0 = *p0*n1n;

    if (n > 2) {
      p1 = d + r;
      p1n = d + r*(n-1);
      m = 2;
      while (p1 < p1n) {
	t=0;
	q = p1 + 1; 
	qn = q + (n-1);
	while (q < qn)
	  t += *q++;

	*(p1 + n-1) = *p1*(m-1) - t;
	p1 += r;
	m++;
      }
    }

    *px++= *d; // m=1
    p1 = d + r;
    p1n = d + r*n;
    while (p1 < p1n) {
      q = p1;
      qn = q + n;
      t = 0;
      while(q < qn) {
	t += *q; 
	q +=1;
      }
      p1 += r;
      *px++ = (double)t;
    }
  }  

  if (d) 
    free(d);
 ERROR:
  return;
}

", "lcg_dist.c")
system("R CMD SHLIB lcg_dist.c")
}

lcg.dist.gen1c <- function(n) {
  d <- c()
  dyn.load(sprintf("lcg_dist%s", .Platform$dynlib.ext))
  for (v in 1:n) {
    r <- .C("lcg_gen_dist1", n=as.integer(v), prob=as.double(rep(0, v)))
    d <- c(d, r$prob[1:v])
  }
  dyn.unload(sprintf("lcg_dist%s", .Platform$dynlib.ext))
  return(d)
}

lcg.dist.gen2c <- function(r) {
  d <- c(rep(0, r*(r+1)/2))
  dyn.load(sprintf("lcg_dist%s", .Platform$dynlib.ext))
  r <- .C("lcg_gen_dist2", n=as.integer(r), prob=as.double(d))
  dyn.unload(sprintf("lcg_dist%s", .Platform$dynlib.ext))
  return(r$prob)
}

# Pr(m; r)
lcg.dist.m.gen <- function(r) {
  d <- c(rep(0, r*(r+1)/2))
  dyn.load(sprintf("lcg_dist%s", .Platform$dynlib.ext))
  r <- .C("lcg_gen_dist_m", n=as.integer(r), prob=as.double(d))
  dyn.unload(sprintf("lcg_dist%s", .Platform$dynlib.ext))
  return(r$prob)
}





